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ABSTRACT 

A uniform approach to construct absorbing artificial boundary conditions for second- 
order linear hyperbolic equations is proposed. The nonlocal boundary condition is given by 
a pseudodifferential operator that annihilates travelling waves. It is obtained through the 
dispersion relation of the differential equation by requiring that the initial-boundary value 
problem admits the wave solutions travelling in one direction only. Local approximation 
of this global boundary condition yields an nth-order differential operator. It is shown 
that the best approximations must be in the canonical forms which can be factorized into 
first-order operators. These boundary conditions are perfectly absorbing for wave packets 
propagating at certain group velocities. 

A hierarchy of absorbing boundary conditions is derived for transonic small pertur- 
bation equations of unsteady flows. These examples illustrate that the absorbing boundary 
conditions are easy to derive, and the effectiveness is demonstrated by the numerical ex- 
periments. 
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1. INTRODUCTION 


Numerical simulations of many physical processes often require calculations of so- 
lutions to partial differential equations on some infinite region. In such problems, it is 
essential to introduce some techniques to restrict calculations to a finite computational 
region. (See [1], [16] and [20].) A traditional method is the coordinates transformation, 
by which an infinite physical region is mapped into a finite computational domain. While 
this technique is effective for some steady state problems, it is inadequate for unsteady 
calculations arising from such applications as seismology, meteorology and transonic fluid 
dynamics. (See [7], [9]). An alternative is to introduce some artificial boundaries to 

obtain a finite region. At these artificial boundaries, some boundary conditions have to 
be prescribed to ensure a unique and well-posed solution. Over the years, there has been 
substantial interest in developing absorbing artificial boundary conditions that eliminate 
the nonphysical reflections, cf. [3], [4], [5], [6], [7], [8], [10], [12], [22], [24] and [25]. 

Using the wave equation as their model, Engquist and Majda [6, 7] proposed a 
pseudo-differential operator as a perfectly absorbing boundary condition. Since the dis- 
persion relation of the wave equation admits two modes, the right-going and the left-going, 
which correspond to the two branches of the square root function ±\/l — s 2 , an equation 
which admits only one of these modes is obtained by choosing a particular branch of the 
square root. This equation is then considered as the dispersion relation of an equation 
containing a pseudodifferential operator. Since a pseudodifferential operator is nonlocal in 
both time and space, this boundary condition is not numerically useful in practice. En- 
gquist and Majda then used Pade-approximation of the square root \/l — s 2 to arrive at 
a hierarchy of local boundary conditions given by differential operators. These boundary 
conditions are perfectly absorbing at normal incidence. 

Higdon [13] proposed a process by which Engquist-Majda boundary conditions can 
be generalized to be perfectly absorbing at some arbitrary angles of incidence. In this 
process, he was able to show that for certain approximations, a differential operator in 
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Engquist-Majda boundary conditions can be put into a canonical form. It can be factorized 
into differential operators of first-order, with each factor annihilating waves at a specific 
angle of incidence. This canonical form has theoretical and practical advantages: it does 
a great deal to simplify stability analysis and numerical implementation [13]. Also this 
characterization of absorbing boundary condition highlights the relation between Engquist- 
Majda boundary conditions and Bayliss-Turkel boundary conditions [5]. 

The purpose of this paper is threefold: 

1) We propose another approach, which is related to that of Engquist and Majda in 
[6] and [7], to construct local absorbing boundary conditions. In their process of rational 
approximation of the square root function, Engquist and Majda have used different strate- 
gies for different equations, namely, the wave equation in [6] and [7] and the transonic 
small disturbance equation in [8]. In our approach, a fairly uniform strategy of rational 
approximation is applied, and the absorbing boundary conditions can be automatically 
generated as long as the dispersion relation is known. 

2) As was pointed out in [13], the factorization theorem of Higdon applies only to 
certain rational approximations for wave equations. We will generalize this factorization 
theorem to show that a larger class of absorbing boundary conditions can be in fact put into 
Higdon’s canonical form. The use of this form will have the advantages both in analysis and 
implementation, especially when higher-order conditions are concerned. When applied to 
such equations as the wave equation and the dispersive Klein-Gordon equation, our method 
yields the same results as established in [13]. 

3) High-order absorbing boundary conditions often involve some optimization pa- 
rameters. Our approach of constructing absorbing boundary conditions provides a natural 
link between these parameters and the group velocities of wave solutions. We will show 
that the optimal absorbing boundary conditions are those which are perfectly absorbing 
at certain group velocities, with each factor annihilating the wave packets propagating 
at a specific group velocity. This physical interpretation is helpful in determining the 
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optimization parameters. 

In section 2, we will derive perfectly absorbing boundary conditions for the travelling 
wave solutions. Its local approximation will be considered in section 3. Also in section 
3, we will describe the general characterization of the absorbing boundary conditions. As 
applications, we will derive a hierarchy of absorbing boundary conditions for the transonic 
small disturbance equations in section 4. Finally, in section 5, we will compare by numerical 
experiments our boundary conditions with Engquist-Majda boundary conditions of [8]. 


2. PERFECTLY ABSORBING BOUNDARY CONDITION 

In this section, we derive a global boundary condition which annihilates travelling 
waves of a general second-order hyperbolic equation. In fact, our global absorbing bound- 
ary condition applies to any linear hyperbolic equations with constant coefficients for which 
a dispersion relation is known. 

We consider the initial-boundary value problem (IB VP): 



o' 

II 

s 

cT 

t > 0, 

xi > 0, 



u(t;0,x_) = g(t; x_), 

t > o, 

x_ e m n_1 , 

(2.1) 


u = 0, 

t < 0, 



where D = 

••> sfc) and x - = (*2,.. 

.,x„). 

P is a homogeneous 

polynomial of 


degree two of n + 1 variables with constant coefficients. Furthermore, we assume P is 
hyperbolic in the direction of t (see [19]) and that g has a compact support. Problem 
(2.1) must be solved over a region which is bounded in x\, e.g., ft = {(®i,. . . ,r n )|ri < a). 

By an absorbing boundary condition Bu = 0 at Xi = a, it is meant that the solution 
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of (2.1) can be well approximated by the solution of the following problem. 


( 2 . 2 ) 


P(—;D)u = 0, t > 0, 0 < Xi < a, 

at 

u(t; 0, x_) = g(t; x_), t > 0, x_ E IR" -1 , 

Bu = 0, = a, t > 0, x_ (E IR" -1 , 

u = 0, t < 0. 

Let (r, £) 6 IR n+1 be the dual of (t,x) and £_ € IR" -1 be the dual of x_. Equation 
(2.1) admits plane wave solutions of the form 


u(<;x) = e <(r<+ * z) , 


( 2 . 3 ) 


for (r, £) 6 IR n+1 , provided the dispersion relation 


= P(ir;iO = 0 , 


( 2 . 4 ) 


is satisfied, where i£ = (i£i , . . . , i£ n ). The group velocity of (2.3) with which the energy of 
the wave packet is propagated is equal to (see [26]) 

C - -Vf • r. 


For a given frequency r € IR and a wavenumber € IR" -1 , we assume that the 
equation in (2.1) admits the plane wave solutions (2.3) which propa,gate in both the positive 
and the negative directions of the xj-axis. Thus for some (r, £_) € IR", (2.4) has two roots 
£j~ 6 IR such that the xj -components of group velocities V *i are positive and negative 
respectively, i.e., 

The travelling waves which propagate in the positive Xj -direction has the form 


ti(<;x) = e i(rt+ tf Zl+ *- , *- ) . 


( 2 . 5 ) 
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It follows that a boundary condition that annihilates the right-going wave solutions 
(2.5) can be written as 

From (2.4), the implicit function theorem leads to 

dr = 

d(i <r T ' 


So the perfectly absorbing boundary condition at x\ = a in the absence of evanescent waves 
is given by 




or 

( 2 . 6 ) 

<7 r a T 

The left hand side of (2.6) can be regarded as the symbol of a pseudo-differential 
operator which can be expressed in terms of a Fourier integral. This operator requires the 
information of the solution away from the boundary xi = a. An alternate definition of the 
boundary condition corresponding to (2.6) can be formulated as follows. 

For u(t-,x ) = define B* by 




(ri-K-x) 


For more general solution, the Fourier transform u(r;ri,£_) of u with respect to t 
and x- can be decomposed into a sum of two terms 


«(r; *,,{_) = / + (r:f-)e‘<. + *' +/"( r;{- )<**'•■. 


(2.7) 


Then the action of B* on u is 


a 


^(r;f+£-)|}/ + (r;«-)e i <' ,+ «f- + «-'->dr<i{. 

^■(r; £f, £-)|}/-(r; f-)e i(r,+{r *>+«-*- >*•*. 


( 2 . 8 ) 
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We can show that if the IBVP (2.1) generates travelling waves only, i.e., the Fourier 
transform g of g satisfies 

supp g C S = {(r, £_) £ IR" | (2.4) has only real roots (2-9) 

then B* defined in (2.8) annihilates the exact solution of (2.1). 

By using the Fourier-Laplace transform, the solution of (2.1) is given by 

u(t;x) = Jj e'( r<+ ^i l ^ l+ ^ - ' x- ^(r; f_)drc?£_. 

In the above integral, = £ + (r; £_) is the root of (2.4) which satisfies Im£^" > 0 for each 
(r;£_) £ (C x IR” -1 with Imr < 0. For r € IR, £^(r;£_) is taken as the limit in the 
lower-half plane Imr < 0. It was shown in [14] that if £ IR, then V^f*) > 0, i.e., the 
positive group velocity in the ri-direction. 

Decompose this solution into (2.7), the result is / + = g and f~ — 0. Then sub- 
stitute it into (2.8). The first integral in (2.8) vanishes by the relation ~^(t; £*,£_) = 
|^(r;£+,f_)| and so does the second integral because f~ = 0. This leads to B*u = 0. 
Thus the solution of (2.1) is a solution of (2.2) if B = B*. 

The solution of (2.2) is also unique when B = B* . The boundary condition B*u = 0 
implies /“(r;f_) = 0, because ^(r;£f,£_) ^ for each ( r :f-)- This 

shows that u is a solution of (2.1) if and only if it is a solution of (2.2). 

Therefore the boundary condition B* = 0, or equivalently (2.6), is a perfectly ab- 
sorbing boundary condition. This boundary condition, however, is nonlocal in both time 
and space due to the presence of the absolute value function in (2.6), hence not useful 
in practice. In order to implement this condition numerically, it must be replaced by 
some local boundary conditions resulting from the rational approximations of the absolute 
value function. Such approximations will be discussed in section 3, but here it is probably 
worthwhile to mention a connection between (2.6) and the perfectly absorbing boundary 
condition derived by Engquist and Majda in [6]. The global boundary condition of [6] 
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involves the square root function which also has to be approximated by some rational 
functions. According to Newman [18], the best rational approximations of the absolute 
value function |x| and the square root y/x have the same order of accuracy. In this sense, 
if (2.6) is approximated by a good rational function, the localized boundary conditions 
resulting from (2.6) are as accurate as those of [6] and [7]. 


3. CHARACTERIZATION OF HIGH-ORDER ABSORBING BOUNDARY 
CONDITIONS 

We now consider the local approximations of the perfectly absorbing boundary con- 
dition (2.6). In designing an absorbing boundary condition two properties must be consid- 
ered: it should minimize the amplitudes of the waves reflected from the artificial boundary 
so that the solution of (2.2) closely approximates the free-space solution of (2.1), and it 
must also be a well-posed condition to guarantee a unique and well-posed solution to the 
differential equation. To study the absorption property of a boundary condition Bu = 0 
with its dispersion relation = 0, we consider the wave solutions of the form 

u(t;x ) = e , ( r ‘ + *i' x i +*-*-> + r e i(rt+(r *!+*-•*-), 

Substitution of u into the boundary condition Bu — 0 at x\ — a yields 

B(r; {+,£_) + r = 0 . 


Solving for r from the above equation, the result is 


= £(r Mt-eri. 

*(r; 

Then the reflection coefficient of the boundary condition Bu = 0 is defined by 


Rb(t', £-) = M = 




B(v, 
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The reflection coefficient is the amplitude of the reflected wave. Therefore a per- 
fectly absorbing boundary condition must yield a reflection coefficient equal to 0 for all 
frequencies, e.g., the boundary condition (2.6). This is, in general, impossible to achieve 
for a local boundary condition. Hence, we hope to build a local boundary condition for 
which Rg is as small as possible. 

A theory to determine the well-posedness of an IBVP has been developed by Kreiss 
[17] and by Sakamoto [21]. The following criterion of well-posedness can be found in [7] 
and [24], 


Well-posedness criterion. 


The IBVP (2.2) is well-posed if and only if 

= 0 , 


(3.1) 


i,£_) = 0, 

has no eigenvalues and generalized eigenvalues. 

An eigenvalue of (3.1 ) is defined as (r; £i,£_) 6 €x (Cx JR n_1 , which satisfies (3.1) 
and Imr < 0 and Im£ i < 0. 

A generalized eigenvalue of (3.1) is defined as (t;£j,£_) 6 IR n+1 with (r; £ 1? £_) / 
(0,0,0), which satisfies (3.1) and 


V Xl (ti) 


< o 

<*v(t;£i,£-) - 


Discussions and interpretations of this criterion can be found in [14] and [24], and 
in [23] for the analogue of difference equations. 

In order to obtain a local boundary condition approximating (2.6), the absolute 
value function in (2.6) must be approximated by a rational function. By making a rational 
transform if necessary, we may assume that V Xl 6 [—1,1]. 

Let r(x) = P g ™l*) °f order (m, n) be a rational approximation to |x| on interval [—1, 1], 
then (2.6) can be approximated by 


°~£i _ r / \ _ Pm(<r^i/<Tr) 

&t q n (o^/o T )' 


(3.2) 
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or 

a^ l af~ 1 q n (— ) - afp m (— ) = 0, d = max{m , n + 1}. (3.3) 

(J T <7 r 

The left hand side of (3.3) is a homogeneous polynomial of degree d in (r, £), hence (3.3) 
is the dispersion relation of a differential equation. 

An obvious choice of r(x ) seems to be the one, r*(®), obtained by Newman in [18], 
which has an accuracy of order e -cv ^. This approximation, however, results in a boundary 
condition that is ill-posed, because r*(0) = 0; hence, (3.3) admits a generalized eigenvalue. 
Although adding a small positive constant to r* would rule out all generalized eigenvalues 
while still maintaining roughly the same accuracy, it is not clear whether the resulting 
boundary condition admits other eigenvalues. Rather than looking for a function approxi- 
mating |® |, we will instead take another approach due to Higdon [13], in which we study the 
necessary form of a well-posed boundary condition with the smallest reflection coefficient. 

Equation (3.2) can be expressed as 



where Q is a polynomial of degree d = max{m, n + 1} with real coefficients, which can be 
factorized as 


d £±4i 


(3.4) 


j= 1 


j=d i+l 


with vj,Wj,Wj G IR. For any polynomial (3.4), there exists a corresponding boundary 
condition whose dispersion relation is given by 


<riQ(—) = 0 . 

& r 


(3.5) 


The boundary condition corresponding to (3.5) will also be referred as boundary 
condition Q of order d. 
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The reflection coefficient of (3.5) is given by 
/? _ TT <T €i(Ci~) ~ 

Q JLwc)-*wer) 

d d ^ 

A (^6 ttt ) - VjQjJJit )) 2 + ) 

x >=1+1 (**(£■) - ^-mg - )) 2 + *&(tn 

with obvious independent variables omitted. 

Now, we assume = ^l(t, £) € [a, /?] for all (r, £_) G «S and £ G IR such that 
cx(r, = 0, and — oo < a < 0 < (3 < +oo; that is, the zj-component of group velocity has 
the lower bound a and the upper bound /?. 

Definition 3.1. Given e : 0 < e < (3, let 


Qd — {Q : real polynomials of degree d}, 

t+ 


5, = {(r,(.) € E" | e < ’f - ) < ,9). 

0>(r; £\f-) 

A boundary condition Q* G is said to be Cf°-best boundary condition of order d 


if 


max f? Q .(r,£_)= min { max J2q(t,£_)}, 
(r,t-)es t QeQi (r,(-)es, 

where Rq is the reflection coefficient of the boundary condition Q. 


£|°- best boundary conditions are the optimal boundary conditions of a fixed degree, 
in the sense that they minimize the reflection coefficient over all travelling waves of prop- 
agating speed in the range [e,/?]. In applications, during the time interval of interest, the 
waves with the slow speed ( V Xl G [0, e)) will not travel far enough to reach the boundary. 

The following definition, due to Wagatha [25], takes into consideration the wave 
modes of every speed in (0, /3]. 


Definition 3.2. 


Let p G C^ 0 p] with p > 0 and 
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An absorbing boundary condition Q* is said to be C l -best boundary condition of 
order d if 


f p(V Xl )R Q -(V Xl )dV Xl = min t p(V Xl )R Q (V Xl )dV Xl . 

Jo QeQi Jo 

The C 2 -best boundary condition is similarly defined by replacing Rq by Rq in the above 
integral. 

The following theorems characterize the best boundary conditions. 


THEOREM 3.3. Let e > 0. If Q* is a well-posed Cf^-best boundary condition, then 

d 

<?•(*) = n<* - ”))• < 3 - 6 > 

j = i 

Furthermore, Vj € [e, /?]. 

Theorem 3.4. Let Q* be a well-posed C 1 ( or C 2 )-best boundary condition. Then (3.6) 
holds and vj 6 (0, /?]. 

The proofs of the theorems 3.3 and 3.4 are identical; we prove Theorem 3.3 only. 


Proof of Theorem 3.3: 

Since Q* is a polynomial of degree d, it can be factorized as 

d 

q*( x ) = n<* - vj ) n k* ~ w i > > 2 + ( 3 - ? ) 

j=i i=d i+i 

(a) We first show that if Q* gives rise to an ££°-best boundary condition, then the 
quadratic factors must disappear in (3.7), i.e., d\ = d . This part of the proof was provided 
by one of the referees. 

We consider a quadratic factor (x — wj) 2 + w 2 and its corresponding reflection coef- 
ficient 

R _ )) 2 + wfcKtf) 

1 (^i (er ) - "MS )) 2 + ) ’ 

= Rjc 2 , 
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where 


■ («+-«,)« +. 8 ? 

J (x - — i«j) 2 + u> 2 ’ 

c = and x* = By assumption, 

a < x~ < 0 < e < x + < 0. 

It is easy to verify, by expanding the both sides, that 


(x + - yjw) + W 2 -) 2 < (X + - Wj ) 2 + Wj, 


and 


This implies 


Rj > 


(x - ywf+ tiff > (x - Wj ) 2 + Wj. 

{x+ - yjw)+wf) 


for all x , 


(x- - ^w) + w))‘ 

The equality in the above holds only if Wj = 0. But because Q* is optimal, the value of 
Rj can not be reduced, hence one concludes xbj = 0. In that case, this quadratic factor is 
in fact a product of two linear factors. This shows that (3.6) holds. 

(b) We now consider a linear factor (x — Vj ) of (3.6) and show that if Vj ^ [e, /?], 
then either the condition (3.6) fails the well-posedness criterion, or it is not optimal. 

If Vj 6 [a, 0], then the condition (3.6) is not well-posed because there exists a gener- 
alized eigenvalue. Now we consider the reflection coefficient corresponding to this factor. 


Rj = 

where x ± and c are defined in part (a). 


X + — Vi 


X — Vi 


C, 


It is easy to verify that if Vj € (— oo,a) U (0, +oo), then 


x 


+ _ 


X — Vi 


> 


x + — 0 


x~ — 0 


for all x , x + . 


Also, if Vj € (0, e), then 


X + 

-Vj 

> 

x + — £ 

x~ 

-Vj 

X~ — £ 


for all x , x + . 


Therefore Vj € [e, 0] for Q* to be optimal. This complete the proof. 


□ 
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Theorem 3.5. Let Q* d be a well-posed C^-best boundary condition. If 

KO;fi + >£-)l ^ kr(r;£~,f_)|, 

for all (r; £_) G S c , then Rq» < 1 for all (t; £_) G S e and 

lim Rq» = 0, for all (r;f_) G S e . 

d—++o o d 

Proof: Let 

d 

Qd(x ) = JJ(x - 20), 

i=i 

then Rqj < 1 and ^ lim Rq a = 0. The theorem follows from Rq» < RQ d . □ 

From Theorem 3.3, the dth order Cf 3 (C 1 or £ 2 )-best condition for the problem (2.1) 
whose dispersion relation is <t(t; £) = 0 must have the form 

U{^D)- V M^D)}u = 0, (3.8) 

j= i 

where D = (^-, . . . , ^|-) , Vj G [e, 0] (or (0, /?)). In particular, the C 1 (or £ 2 )-best bound- 
ary conditions can be generated automatically with the help of a symbolic manipulation 
program in the following manner. First, the derivatives cr ^ , o r are calculated from the 
dispersion relation <r(r; £) = 0. Then the parameters vj can be approximated through a 
numerical process to minimize the integral in Definition 3.2. Finally, analysis for well- 
posedness can be carried out for each factor of (3.8), which is simple as shown in the next 
section. 

4. ABSORBING BOUNDARY CONDITIONS FOR A CLASS OF TRAN- 
SONIC SMALL DISTURBANCE EQUATIONS OF UNSTEADY FLOWS 

In their paper [8], Engquist and Majda derived absorbing boundary conditions up 
to the second-order for the unsteady transonic small disturbance equation of low reduced 
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frequencies. Their derivation, based on the framework of their earlier work [7], did not, 
however, provide a uniform approach in constructing the local boundary conditions. Dif- 
ferent strategies had to be used to approximate the square root while handling the sidewall 
conditions. Furthermore, if a higher-order approximation is required, their method leads 
to two unpleasant difficulties: the determination of well-posedness and the selection of the 
parameters. Although general guidelines for selecting the parameters were given, these 
guidelines do not have a direct physical interpretation. 

In this section, we will recreate the artificial boundary conditions for the transonic 
small disturbance equation of low reduced frequencies treated in [8]. From our approach, 
a connection between absorption properties and group velocities of the interior distur- 
bance is revealed. Boundary conditions for the transonic small disturbance equation of full 
frequencies will also be considered. 

4.1 Equation of Low Reduced Frequencies 

By a standard frozen coefficient theory, we may assume that the equation is linear 
with constant coefficients and the flow is subsonic, which is the case in farfield. Thus we 
consider 


2<f> xt = K*<f> xx + <f> yy , K*> 0 (4.1) 

1*1 < ll/I < b. 

Let (r, £, rj) be the dual of ( t,x,y ), then the dispersion relation of (4.1) becomes 

a(r,(,r t ) = 2T(-K*t 2 -T 1 2 =0. (4.2) 

We first consider the sidewalls. 


Sidewall condition 

The y- component of the group velocity of the plane wave e'( rt +? I + T ?J') j s equal to 

T/ _ _ 2*7 _ rj 

y — — — — — - . 


2 £ 


£ 
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Dividing through (4.2) by £ 2 , we have 

(|) 2 = 2 j-K'. (4.3) 

From this relation, we are able to determine the range of group velocity. Since the right- 
hand side can be made arbitrarily large, the group velocity is unbounded, i.e., 

V y e (-oo,+oo). 

Even though the propagating speed in the y-direction may be infinity, it can be 
justified to consider only the group velocities in a bounded interval, for the following 
reasons. 

The propagating speed of a plane wave in the y-direction becomes infinity only if 
the wave is parallel to z-axis, i.e., f = 0. In this case, as is easy to verify, the propagating 
speed downstream is also infinite and is much faster than the y-component because 

^jr — ► 0 as £ — ► 0. 

V x 

This implies that the disturbance with infinite y-direction speed will actually strike the 
downstream boundary first. So in designing a sidewall boundary condition, we need only 
to consider those waves with finite propagating speed in the y-direction. Thus, we may 
assume there is a /? > 0 such that 

v.eHM. 

For a given (r, f) € IR 2 , equation (4.2) has two roots tj + = —T]~, corresponding to 
the two travelling wave modes, with 

o < v»(-? + ) = -non- 

Thus as a direct consequence of Theorem 3.3, an ££°-best boundary condition of order d 

at sidewall y = b (resp., y = — b ) must have the form 

d 

Qd(x) = JJ(x - Vj), (4.4) 

>= i 
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with vj E [e, /3] (resp., Vj E [— /3, — e]), e > 0. If the reflection coefficient of (4.4) is denoted 
by Rq, then it follows from Theorem 3.5 that 

- r q( t iO< 1 for all (r,f) € S e , 


where S e is defined in the Definition 3.1. 

The theory of section 3 does not assure that the boundary condition given by (4.4) 
is well-posed; it has to be determined by the criterion given at the beginning of section 3. 
In order to prove the well-posedness, we need only to show that each factor in (4.4) yields 
a well-posed boundary condition. Thus, let us consider the factor ( x — Vj). The dispersion 
relation of the corresponding boundary condition is 



or 


7] + Vj C = 0. (4.5) 

Let us consider the mutual solutions of (4.2) and (4.5). For any £ E IR, the solution 
r) of (4.5) is always real, thus an eigenvalue does not exist. Also, for any £ E IR, (£, r/ + ) 
will never satisfy (4.5) since vj < 0. This implies there is also no generalized eigenvalues, 
so the boundary condition is well-posed. 

From this, the best boundary condition at the sidewalls should have the form 

n<& + "5>'-°- < 4 - 6 ) 

with vj < 0 for y = — b and vj > 0 for y = b. The parameters Vj's must be tuned to make 
the reflection coefficient small. Details in the determination of parameters will be given in 
section 5. In general, if we know a 'priori the main interior disturbances propagate with 
certain speeds in the y-direction, these disturbances can be exactly annihilated by tuning 
Vj to those speeds. 
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In [8], = —y/K* was chosen for the first-order condition, which annihilates exactly 

the waves travelling in the y-direction only, i.e., with zero speed in the ar-direction. 

Upstream condition 

The ar-direction of group velocity is given by 

V, - a = _jr + ' 


By employing the previous techniques, we find that 




It follows from Theorem 3.3 that the ££°-best boundary condition of order d at 
x = —a has the form (which is well-posed): 


n (!-<*• +*>£>-«• 


(4.7) 


i= i 


where vj 6 [— — e] . Its reflection coefficient is smaller than 1 by Theorem 3.5. For the 

first-order condition, v\ = was chosen in [8], the corresponding boundary condition 

annihilating waves travelling at the fastest speed upstream. 

The second-order condition of [8] is given by 


$ xt j{* < ^ >tt 2^ yy ~ 


which is also a special case of (4.7). In fact, if one solves for <f> yy from (4.1) and substitutes 
it into the above equation, the result is 

2 . K * 

*<Pxt Ytt g tZX — 0* 

K* 


This reduces to (4.7) with v\ = V 2 = — 
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Downstream condition 

At the downstream x = a, the ££°-best boundary condition also has the form (4.7) 
but with Vj > 0. Condition (4.7) is well-posed at x = a if Vj > 0, and its reflection 
coefficient is equal to 


TH \ — TT ( Yl ~ — Ylh/1 ^*^r) 2 ^ 

T ’ “ JJ U + (K- + OiW-K'Wj ' 

It is clear that R(r,rj) < 1 if yjl — K*^) 2 ^ 0, i.e., V x (£) ^ 0. 

Since the interior disturbance might contain certain modes which have infinite down- 
stream speed, one of the parameters can be tuned to annihilate these wave modes, e.g., 
vi = +00 in the first-order condition. In this case, the reflection coefficient becomes 

Rl{ ,n) ~ l + Jl-K'W 


and the corresponding boundary condition is obtained by taking the limit v\ — * -f-oo: 


— <j> = 0, at x = a. 
ox 


This is exactly the boundary condition given in [8]. For higher order, we propose 


a — 1 

n( 


i d 
K* + vj dt 


d_ 

dx 



= 0 . 


(4.8) 


4.2 Equation of Full Frequencies 

We now consider the full frequency transonic small disturbance equation 

(f>tt d" 2 (f> x t = A <f) xx d - fiyy- 

Its dispersion relation is 

v(r , £, 77) = r 2 d- 2£r - 2 - rj 2 = 0. (4.9) 


Sidewall condition 
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The y-component of group velocity is equal to 



*7 

t + C 


It can be shown that 


v, € [- 1 , 1 ], 


and the two travelling wave modes satisfy 


77 + = -t] and |V r y (^ + )| = \V y (rj )|, |<r r (T 7 +)| = |oy(t 7 )|. 


Applying Theorem 3.3 and Theorem 3.5, the ££°-best condition has the form 

d 

q*(x ) = n (* - 

i = i 

Vj G [—1, —e] for y = —b and vj G [e, 1] for y — b. Furthermore Rq* < 1. Therefore, the 
sidewall conditions may be chosen as 


P( 

i = 1 



d_ 

dt 



o. 


The condition is well-posed. 


Upstream and downstream conditions 

From the dispersion relation (4.9), we find 

V x = 1 — (K* + 1 ) — -j-t, 

T + ( 


and 


v x 6 [i- Vi + k*, 1 + Vi + K*}. 


It follows from Theorem 3.3 that the absorbing boundary conditions can be taken as 


n ( (i - ">4 ■ (A " + vi) h )4, = °’ <4io) 
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where vj G [1 — \/l + A'*, — e] for upstream condition and Vj € [«:, 1 + y/l + A*] for down- 
stream condition. 

It can be verified directly that the reflection coefficients associated with (4.10) for 
both upstream and downstream are strictly less than 1, except for the mode with zero 
x-component group velocity. 

5. NUMERICAL COMPUTATIONS 

In this section, we discuss the numerical implementation and present some compu- 
tational results involving the absorbing boundary conditions derived in the last section. 
For the purpose of illustrating the efficiency of these conditions, we compare the com- 
putational results obtained by using the boundary conditions of this paper with those 
developed by Engquist and Majda in [8], both applied to the linear equation (4.1) in the 
region 17 = [—1,1] x [—1,1]. We note that the first-order conditions derived from the last 
section are the same as those derived in [8]. 

Specifically, we compare the following sets of conditions: 

Upstream: 

First-order: 

4>t + v<j> x — 0. 

Second-order: 


. 2 1 
Yxt — <PU + 2 < es/y — 0, 

(E k M) 


(5.1) 


Downstream: 

First-order: 

4> x — 0 . 
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Second-order: 


<t>xx — Ol 


f— - 

\ A'* + v dt 


d_ 

dx 



= 0 . 


Sidewall ( y — 1): 
First-order: 


Second-order: 


<f> y + v<f> x = 0. 


<f>ty + r 2 <t>xy + ri<f>xt — 0 , 

/ d d \ ( d 

ya-y +v '-^)\d^ +V2 dir = 0 - 


(E & M) 
(5.2) 


(E & M) 
(5.3) 


5.1 Parameters 

Theoretically, the parameters in the above boundary conditions should be determined 
by the minimization properties stated in the definitions 3.1 and 3.2. Such process, however, 
is hardly necessary in practice, nor will it always yield the best result for all problems, 
because the definition of the or £ 2 -best conditions itself depends on the arbitrary 
choice of e or the weight p. In the present study, we use the following procedure in which 
the parameters can be easily determined. 

For the upstream conditions, the parameters must satisfy v € [-^-,0). The first- 
order condition is tuned to the waves travelling most rapidly upstream by letting v = — 

In the second-order condition (5.1), the two parameters are chosen as «i = — and 
i >2 = — the mid-point of [-^-,0). 

The parameters in the downstream conditions must be in (0, +oo), which is un- 
bounded. However, a transformation exists to map (0, +oo) to the bounded interval (0, 1), 
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hence the similar method as above can be used. Recall from section 4 that a factor of 
reflection coefficient of the downstream condition is 


Rj = 


Vj - (K- + VjWl - 
VJ + ( K * + Vj)y/i -*•(?)*’ 

- V' - K-ar 
~ X^ + V^T^W' 


whore S (°. !)• 

Thus the first-order downstream condition is obtained by letting = 1, i.e., 

v = +oo. The remaining parameter v in the second-order condition (5.2) is determined by 
requiring K ? +v = |, therefore, v = K*. 

In the case of sidewall, we have to deal with the genuine unbounded interval. For 
the reasons mentioned in section 4, v = y/K* in the first-order sidewall (y = 1) condition. 
Parameters tq, v 2 in (5.3) are determined by minimizing the integral 



e -«(*-vW R 2 ( x y dx , 


where, by section 4, R 2 = (f+^)(f+g)- 

By a change of variables, the above integral becomes 



g — 6(x— 1)* / ® h) 2 ( X ~ d2 ) 2 dx 

' x + Vi a; + v 2 ' 


where 0 = Vk?' a = aK* and Vj = Vk '■ v\ and v 2 can be computed by a numerical 
integration and optimization algorithm. For 0 = 10. and a = 0.1, we found, correct to 
the first digit, v\ = 3.0 and v 2 = 0.4, hence tq = 3.0 y/K* and i>;» = 0Ay/K* in (5. 3). This 
particular choice of parameters yields a reflection coefficient R 2 < 0.05 for all wave modes 
with the group velocities in the range [0.25, 5.0] (see [15]). 

Finally, in this study, r x = 6.4641 y/K* and r 2 = 5.4641RT* in (E & M) second-order 
condition, according to one of the guidelines of [8]. 
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5.2 Discretization 

Equation (4.1) is discretized by a semi-implicit scheme described in [2], where the 
ar-direction is treated explicitly and the y-direction is kept implicit. The resulting differ- 
ence equation are solved by marching both in the x- and ^-direction. At a given time 
level, tridiagonal equations are solved along the y-directions, successively marching in the 
downstream direction. For stability, the C.F.L. condition K* ^ < 2 must be satisfied. 
The scheme is nondissipative when K* ^ = 2. 

For the first-order boundary conditions and the second-order Engquist-Majda (E & 
M) boundary conditions, we use exactly the discretization suggested in [8]. Since each 
factor in (5.1), and (5.3) is in the form of the first-order operator of the corresponding 
first-order boundary condition, the discretizations of (5.1) and (5.3) are obtained by apply- 
ing twice the difference operators which discretize the corresponding first-order boundary 
conditions of Engquist and Majda. For example, in [8], the boundary condition 


(f)y -f“ V(j)^ 0, 


is discretized by 


I - K 
Ay 


-x r r-i 


(/+*->)}*;£, =o, 


where J, K are shift operators in x and y respectively, = <j>J +1 k , K<f>J k = 4>j k+1 - 

By applying this twice, the difference equation for (5.3) becomes 

* {^— 1 + • / ' 1 ) + "’^^( / + *~‘)Kth = 0 ' (5-4) 

Condition (5.2) is discretized by an explicit scheme, backward in x. 

The stability of these difference schemes can be determined by the GKS theory [11], 
or by numerical experiments. Because of the factorization form, the difference schemes 
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for the higer-order boundary conditions are stable if and only if, by the GKS theory, the 
schemes for the corresponding first-order conditions are stable. The discretization of the 
first-order boundary conditions is discussed in detail in [8], and although no analysis has 
been attempted, these schemes are found to be stable by numerical experiments, even for 
the nonlinear problem. 

Since the interior equation is solved line by line marching in the downstream direction, 
the boundary condition (5.4) which requires values of <f> at j — 1, j — 2 cannot be applied to 
the immediate neighboring point of the corner of the upstream and the sidewall boundaries. 
This difficulty is overcome by using the following condition at the upstream corner: 


a{x) + ° J; V + (1 - a(i)) (| + ”■ £) (I + 1,2 £) ■ * = °- 


where 0 < a < 1. a decreases smoothly from a = 1 at the upstream boundary to 0. In 
our study, the transition region of a from 1 to 0 extends to 2 grid points. 

There is no difficulty at the downstream comers, since these corner points are not 
involved in the computation. 


5.3 Energy 

In the following computations, we will use the energy of solutions as a means of 
measurement. The energy 


E(n) = £**{ thhtil ± +J±hLl tlm - +h } 2 AlAy (5.5) 

j, k 


+ £ { + - * lt } 2 AxAy, 

h k y 


is the discrete version of 


m = ffr ||1 (a'* <t>l + <f>l)dxdy. 

It is easy to verify that if <t> satisfies (4.1) in |x| < 1 and \y\ < 1, and the perfectly reflecting 
boundary conditions 

<t> = 0, at x — ±1 and y = ±1, (5.6) 
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then the energy of solutions is conserved, that is, 

J 

— E(t) = 0, for all t. (5.7) 

at 

In order to examine the effects of each absorbing boundary condition, without be- 
ing influenced by the presence of the other boundary conditions, we will use a strategy 
where the absorbing boundary condition will be prescribed at one boundary only in each 
experiment, and the perfectly reflecting boundary condition (5.6) will be imposed on the 
remaining boundaries. In the first group of calculations, the energy (5.5) of solutions in 
the domain Q, = [—1, 1] x [—1, 1] will be computed. Therefore, the rate of decrease of the 
energy will demonstrate the ability of the test boundary condition to radiate the energy 
away from the computational domain, because by (5.7), energy of the solution will be 
totally reflected back from the other three perfectly reflecting boundaries. 

In all computations, K * = 1 is used in the equation (4.1), and the region Q = 
[— 1, 1] x [—1,1] is covered by a 60 x 40 uniform mesh, with a CEL condition of 1.997. 
For an initial value, we use a pulse of compact support, a piece of radially symmetrical 
sine function. Energy of solutions is calculated after every 10 time steps, with a total 
computation of 230 time steps. The initial pulse will soon spread out in a parabolic 
wavefront [8], and total reflections will occur at the three perfectly reflecting boundaries, 
so after some time, the disturbance that strikes the test boundary will consist of wave 
packets with a fairly large spectrum of frequencies and wave numbers. 

Figure 1 shows the energy E°(n), given by (5.5), of the solution of (4.1) with the 
perfectly reflecting boundary conditions (5.6) at different time steps. The decrease in the 
energy is due to the dissipation of the difference scheme. In order to eliminate the effect 
of this dissipation in our study, the energy of the solutions for other calculations will all 
be scaled by E°(n), i.e., the energy 

c(n) = ( J£o) * 100 ' 
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will be used for other calculations. Therefore e(n) = 100 if the perfectly reflecting boundary 
condition <f> = 0 is prescribed on all boundaries. 

The results of calculations involving boundary conditions on the upstream, down- 
stream and sidewall (y = 1) boundaries are presented by graphs in Fig. 2, Fig. 3 and Fig. 4, 
respectively. The curves of the energy e(n) all follow one pattern: during the short period 
of the initial time (about 30 time steps), the three (one first-order and two second-order) 
boundary conditions produce roughly the same results, but they differ significantly after 
a long time period (after about 130 time steps). 

In the second group of calculations, we measure the energy of the reflected waves 
of different absorbing boundary conditions. The reflected wave is obtained by comparing 
the solution of an absorbing boundary condition with the free-space solution which is 
calculated in a larger computational area. For example, in the computations involving 
upstream boundary conditions, the free-space solution is calculated in the region U fi 
as illustrated in Fig. 5. Then an absorbing boundary condition is used on the boundary 
x = — 1, and the corresponding solution is computed in the region Q. The difference of these 
two solutions in 0 can be considered as the reflection caused by the absorbing boundary 
condition. Similar methods are used for the downstream and the sidewall boundaries. 

The results for the energy of reflected waves are shown in Fig. 6, Fig. 7 and Fig. 8. 
As before, the energy is calculated at every 10 time steps. In each graph, the curves are 
scaled so that the maximum energy in the first-order boundary condition is 100. 

These numerical results clearly show the improved performance of the second-order 
absorbing boundary conditions over the first-order conditions. In long time computations, 
the second-order conditions (5.1), (5.2) and (5.3) produce less reflections than the corre- 
sponding (E & M) conditions. However, the second-order (E fc M) conditions generate 
smaller reflections during the very short period of the initial time because they are tuned 
to the waves travelling most rapidly. These observations agree with the analysis of the 
reflection coefficients given in the section 4. 
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Fig. 1: Energy of the solution, with the perfectly reflecting condition on all four sides. 
The decrease in the energy is due to the dissipation of the scheme. 

Fig. 2: Energy of the solutions: Upstream boundary, vi = — v 2 — — in condi- 
tion (5.1). 

Fig. 3: Energy of the solutions: Downstream boundary, v — K* in condition ( 5 . 2 ). 

Fig. 4: Energy of the solutions: Sidewall (y = l) boundary. r x = 6.4641 y/K', r 2 
5.4641 K* in (E & M), and = 3.0 \/K*, v 2 = 0 A\fK* in condition (5.3). 
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Fig. 5: Large computational area for free-space solution calculation: Upstream bound- 
ary. 

Fig. 6: Energy of the reflected waves: Upstream boundary. 

Fig. 7: Energy of the reflected waves: Downstream boundary. 

Fig. 8: Energy of the reflected waves: Sidewall (y = l) boundary. 
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